{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numba import jit\n",
    "\n",
    "@jit(nopython=True)\n",
    "def deltaE(array, x, y):\n",
    "    deltaE = 2 * array[x, y] * (\n",
    "        array[(x - 2) % (2*L), y] * array[(x - 1) % (2*L), (y - 1) % L] +\n",
    "        array[(x + 2) % (2*L), y] * array[(x + 1) % (2*L), (y - 1) % L] +\n",
    "        array[(x - 1) % (2*L), (y + 1) % L] * array[(x + 1) % (2*L), (y + 1) % L]\n",
    "    )\n",
    "    return deltaE\n",
    "\n",
    "@jit(nopython=True)\n",
    "def w(deltaE, T):\n",
    "    if deltaE <= 0:\n",
    "        return 1.0\n",
    "    else:\n",
    "        return np.exp(-deltaE / T)\n",
    "\n",
    "@jit(nopython=True)\n",
    "def E_array(array):\n",
    "    E = 0\n",
    "    for i in range(2*L):\n",
    "        for j in range(L):\n",
    "            E -= array[i, j] * array[(i - 2) % (2*L), j] * array[(i - 1) % (2*L), (j - 1) % L]\n",
    "    return E\n",
    "\n",
    "@jit(nopython=True)\n",
    "def Metropolis(array, T):\n",
    "    x = np.random.randint(0, 2*L)\n",
    "    y = np.random.randint(0, L)\n",
    "    while array[x, y] == 0:\n",
    "        x = np.random.randint(0, 2*L)\n",
    "        y = np.random.randint(0, L)\n",
    "    delE = deltaE(array, x , y)\n",
    "    if np.random.rand() < w(delE, T):\n",
    "        array[x, y] = -array[x, y]\n",
    "    return\n",
    "\n",
    "def main():\n",
    "    for T in np.arange(Tmin, Tmax, deltaT):\n",
    "        array = -np.ones((2*L, L), dtype=np.int8)\n",
    "        #array = np.random.choice([-1, 1], size=(2*L, L))\n",
    "        array[1::2, 1::2] = 0\n",
    "        array[::2, ::2] = 0\n",
    "        E_temp = np.zeros(N1)\n",
    "        M_temp = np.zeros(N1)\n",
    "        for i in range(N1):\n",
    "            Metropolis(array, T)\n",
    "            E_temp[i] = E_array(array)/(L*L)\n",
    "            M_temp[i] = 2*np.abs(np.mean(array))\n",
    "        E.append([T, E_temp])\n",
    "        M.append([T, M_temp])\n",
    "    return\n",
    "\n",
    "L = 40\n",
    "N1 = 10000 * L * L\n",
    "Tmin = 0.2\n",
    "Tmax = 1.4\n",
    "deltaT = 0.2\n",
    "\n",
    "E = []\n",
    "M = []\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    main()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Extract the energy values for each temperature\n",
    "energies = [item[1] for item in E]\n",
    "temperatures = [item[0] for item in E]\n",
    "\n",
    "# Plot the energy values for all temperatures\n",
    "for i in range(len(temperatures)):\n",
    "    x = np.arange(len(energies[i]))\n",
    "    #x = temperatures[i] * np.log(x)  # Multiply x-axis by T * log(t)\n",
    "    x = np.log(x)\n",
    "    plt.plot(x, energies[i], label=f'Temperature: {temperatures[i]}')\n",
    "\n",
    "plt.xlabel('log(t)')\n",
    "plt.ylabel('Energy')\n",
    "plt.title('Energy vs T * log(t)')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
